/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		P52.cpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-20 04:15
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-09-03 19:30
*
*   Description:
*
=============================================================================*/
#include "../../include/fem.hpp"
#include "../../include/implicit_solver.hpp"
#include "../../include/incremental_load.hpp"
#include "../../include/model/model.hpp"
#include "../../include/solid.hpp"
#include "../../include/static_equation.hpp"

LOGGING_REGISTER()

int main(int argc, char* argv[]) {

  using namespace geoxfem;

  using MethodPolicy   = model::fem;
  using CouplePolicy   = model::solid;
  using DimensionType  = geometry::plane_strain;
  using EquationPolicy = model::static_equation;
  using SolvePolicy    = model::implicit_solver;
  using LoadPolicy     = model::incremental_load;
  using Model =
    geoxfem::model::Model<MethodPolicy, CouplePolicy, DimensionType,
                          EquationPolicy, SolvePolicy, LoadPolicy, mpl1>;

  std::string input_file_name = "test/test_P52/P52.json";

  /*  detail, trace, process, info, warning, debug, error }; */
  geoxfem::log::log_init(geoxfem::log::severity_level::error);

  io::DataParser<io::json> parser;
  parser.read(input_file_name);

  /* io::DataExtractor<Model> extractor; */
  /* extractor.init(parser); */

  io::DataWriter<Model, io::vtu> writer;
  writer.init(parser.getPath(), parser.getFileName());

  io::DataPloter ploter;

  Model model;
  /* model.append(extractor); */
  model.append(&writer);

  model.init(parser);
  model.buildProblem();
  model.solveProblem();

  for(auto& v : shapesIn(*model.getMesh(), 0))
    if(v->ifBoundary())
      v->plot(ploter, " w p pt 2 lc rgb \"blue\"");
  ploter.appendSetting("set view 0,0");

  MatrixXR correct_ans(12, 2), ans(12, 2);
  correct_ans << 0.0000E+00, -0.8601E-05,  //
    0.0000E+00, -0.4056E-05,               //
    0.0000E+00, 0.0000E+00,                //
    -0.5472E-07, -0.3771E-05,              //
    0.1225E-05, -0.1910E-05,               //
    0.0000E+00, 0.0000E+00,                //
    0.1792E-06, 0.5860E-06,                //
    0.6280E-06, 0.1707E-06,                //
    0.0000E+00, 0.0000E+00,                //
    0.0000E+00, 0.1129E-06,                //
    0.0000E+00, 0.1066E-06,                //
    0.0000E+00, 0.0000E+00;

  for(auto& v : shapesIn(*model.getMesh(), 0))
    ans.row(v->getIndex()) << v->toVertex().getDof("x0")->getCurrValue(),  //
      v->toVertex().getDof("x1")->getCurrValue();

  std::vector<std::tuple<_real, _real, _real, _real, _real, _real>> data;
  for(auto vertex : shapesIn(*model.getMesh(), 0)) {
    auto  gp    = vertex->toVertex().getGauss();
    _real x     = gp.x();
    _real y     = gp.y();
    _real value = gp.getCurrSigma()(0, 0);
    data.emplace_back(std::make_tuple(x, y, value, 0, 0, 0));
    /* plt.plotPoint(x, y, value); */
  }
  for(auto shp : shapesIn(*model.getMesh(), 2))
    for(auto& gp : gpsIn(*shp)) {
      _real x     = gp.x();
      _real y     = gp.y();
      _real value = gp.getCurrSigma()(0, 0);
      data.emplace_back(std::make_tuple(x, y, value, 0, 0, 0));
      /* plt.plotPoint(x, y, value); */
    }

  /* plt.addSetting("set pm3d at b"); */
  /* plt.addSetting("set dgrid3d 50,50,50"); */
  /* plt.addPlotCmd("\'-\' using 1:2:3 w l"); */
  /* plt.addData(data); */
  /* plt.print(); */
  /* plt.show(); */
  /* model.getMesh()->plot(plt); */
  /* plt.plotContour(data); */

  arithmetic::Stress<DimensionType> shp0gp1, correct_shp0gp1,  //
    shp0gp2, correct_shp0gp2,                                  //
    shp0gp3, correct_shp0gp3,                                  //
    shp0gp4, correct_shp0gp4;

  /* auto shp0 = model.getMesh()->getShape(2, 0); */
  /* shp0gp3 = shp0->getGaussPoint(3).getCurrSigma(); */
  correct_shp0gp3 << -0.4299E+00, 0.1431E+00, 0,  //
    0.1431E+00, -0.1058E+01, 0,                   //
    0, 0, 0;
  /* shp0gp3 = shp0->getGaussPoint(3).getCurrSigma(); */
  correct_shp0gp2 << -0.2511E+00, 0.8629E-01, 0,  //
    0.8629E-01, -0.6411E+00, 0,                   //
    0, 0, 0;
  /* shp0gp3 = shp0->getGaussPoint(3).getCurrSigma(); */
  correct_shp0gp4 << -0.3304E+00, 0.8352E-01, 0,  //
    0.8352E-01, -0.1016E+01, 0,                   //
    0, 0, 0;
  /* shp0gp3 = shp0->getGaussPoint(3).getCurrSigma(); */
  correct_shp0gp1 << -0.1516E+00, 0.2667E-01, 0,  //
    0.2667E-01, -0.5985E+00, 0,                   //
    0, 0, 0;

  arithmetic::Stress<DimensionType> shp1gp1, correct_shp1gp1,  //
    shp1gp2, correct_shp1gp2,                                  //
    shp1gp3, correct_shp1gp3,                                  //
    shp1gp4, correct_shp1gp4;

  /* auto shp1 = model.getMesh()->getShape(2, 1); */
  correct_shp1gp3 << -0.2856E+00, 0.8499E-01, 0,  //
    0.8499E-01, -0.9141E+00, 0,                   //
    0, 0, 0;
  correct_shp1gp2 << -0.1427E+00, 0.1394E+00, 0,  //
    0.1394E+00, -0.5806E+00, 0,                   //
    0, 0, 0;
  correct_shp1gp4 << -0.3808E+00, 0.3735E-01, 0,  //
    0.3735E-01, -0.9549E+00, 0,                   //
    0, 0, 0;
  correct_shp1gp1 << -0.2379E+00, 0.9177E-01, 0,  //
    0.9177E-01, -0.6215E+00, 0,                   //
    0, 0, 0;

  arithmetic::Stress<DimensionType> shp2gp1, correct_shp2gp1,  //
    shp2gp2, correct_shp2gp2,                                  //
    shp2gp3, correct_shp2gp3,                                  //
    shp2gp4, correct_shp2gp4;

  /* auto shp2 = model.getMesh()->getShape(2, 2); */
  correct_shp2gp3 << -0.1513E+00, 0.6412E-01, 0,  //
    0.6412E-01, -0.3680E+00, 0,                   //
    0, 0, 0;
  correct_shp2gp2 << 0.2765E-03, 0.1010E+00, 0,  //
    0.1010E+00, -0.1429E-01, 0,                  //
    0, 0, 0;
  correct_shp2gp4 << -0.2159E+00, 0.1359E-01, 0,  //
    0.1359E-01, -0.3957E+00, 0,                   //
    0, 0, 0;
  correct_shp2gp1 << -0.6432E-01, 0.5050E-01, 0,  //
    0.5050E-01, -0.4197E-01, 0,                   //
    0, 0, 0;

  arithmetic::Stress<DimensionType> shp3gp1, correct_shp3gp1,  //
    shp3gp2, correct_shp3gp2,                                  //
    shp3gp3, correct_shp3gp3,                                  //
    shp3gp4, correct_shp3gp4;

  /* auto shp3 = model.getMesh()->getShape(2, 3); */
  correct_shp3gp3 << -0.2331E+00, 0.1477E+00, 0,  //
    0.1477E+00, -0.4231E+00, 0,                   //
    0, 0, 0;
  correct_shp3gp2 << -0.9447E-01, 0.1211E+00, 0,  //
    0.1211E+00, -0.9964E-01, 0,                   //
    0, 0, 0;
  correct_shp3gp4 << -0.1867E+00, 0.1015E+00, 0,  //
    0.1015E+00, -0.4032E+00, 0,                   //
    0, 0, 0;
  correct_shp3gp1 << -0.4805E-01, 0.7493E-01, 0,  //
    0.7493E-01, -0.7974E-01, 0,                   //
    0, 0, 0;

  arithmetic::Stress<DimensionType> shp4gp1, correct_shp4gp1,  //
    shp4gp2, correct_shp4gp2,                                  //
    shp4gp3, correct_shp4gp3,                                  //
    shp4gp4, correct_shp4gp4;

  /* auto shp4 = model.getMesh()->getShape(2, 4); */
  correct_shp4gp3 << 0.1061E-02, -0.4210E-01, -0,  //
    -0.4210E-01, 0.7274E-01, 0,                    //
    0, 0, 0;
  correct_shp4gp2 << -0.2619E-01, -0.2217E-01, -0,  //
    -0.2217E-01, 0.9144E-02, 0,                     //
    0, 0, 0;
  correct_shp4gp4 << -0.3382E-01, -0.3302E-01, -0,  //
    -0.3302E-01, 0.5779E-01, 0,                     //
    0, 0, 0;
  correct_shp4gp1 << -0.6108E-01, -0.1308E-01, -0,  //
    -0.1308E-01, -0.5806E-02, 0,                    //
    0, 0, 0;

  arithmetic::Stress<DimensionType> shp5gp1, correct_shp5gp1,  //
    shp5gp2, correct_shp5gp2,                                  //
    shp5gp3, correct_shp5gp3,                                  //
    shp5gp4, correct_shp5gp4;

  /* auto shp5 = model.getMesh()->getShape(2, 5); */
  correct_shp5gp3 << -0.4855E-01, 0.3616E-01, 0,  //
    0.3616E-01, 0.1373E-01, 0,                    //
    0, 0, 0;
  correct_shp5gp2 << -0.5281E-01, 0.8267E-02, 0,  //
    0.8267E-02, 0.3775E-02, 0,                    //
    0, 0, 0;
  correct_shp5gp4 << 0.2640E-03, 0.3758E-01, 0,  //
    0.3758E-01, 0.3465E-01, 0,                   //
    0, 0, 0;
  correct_shp5gp1 << -0.4001E-02, 0.9688E-02, 0,  //
    0.9688E-02, 0.2469E-01, 0,                    //
    0, 0, 0;

  return 0;
}
